%% ��������,��load,k,cixu,save
load data80_deal;
load data80;
clear arbc qx_narx;
k=[0.92927,0.049232];
ri=122;
cixu=4;%Ҫ�ڼ��õ�,��Ϊͬһ�߶�,ͼ�����,�����Ҳ�ͬ���
%% ͨ��
%% ���Ȳ���AR(1)ģ��
qxs=reshape(qx,96,ri);
cfs=reshape(cf,96,ri);
arbc(1,:)=es(end,1:ri-1)*k(1)+k(2);
for i=2:96
    arbc(i,:)=arbc(i-1,:)*k(1)+k(2);
end
qx_ar=qxs(:,2:end)+arbc;
qx_ar(qx_ar<0)=0;
error_ar=cfs(:,2:end)-qx_ar;
for i=1:96
    e_ar(i,:)=mae(error_ar(i,:));
end

%% ���Է������Իع�ģ��
for i=1:ri-1
    phInputs_c=num2cell(qx(96*i:96*(i+1))');
    PhTargets_c=num2cell(e(96*i:96*(i+1))');

    [p1,Pi1,Ai1,t1] = preparets(narx_net_closed,phInputs_c,{},PhTargets_c);
    % �������
    yp1 = narx_net_closed(p1,Pi1,Ai1);
    qx_narx(:,i)=qxs(:,i+1)+cell2mat(yp1)';
%     figure;
%     plot([cell2mat(yp1)' cell2mat(t1)'])
end
qx_narx(qx_narx<0)=0;
error_narx=cfs(:,2:end)-qx_narx;
for i=1:96
    e_narx(i,:)=mae(error_narx(i,:));
end

figure;
plot(e_origin,'-*','Markersize',6);
hold on
plot(e_ar,'r-+','Markersize',6);
hold on
plot(e_narx,'m-o','Markersize',6);
% hold on
% plot(e_narxm,'k-.','Markersize',6);
set(gca, 'FontSize', 30);
% hl=legend('Error-Origin','Error-AR(1)','Error-Narx');
% set(hl,'Orientation','horizon');
% set(hl,'Box','off');
ylabel('MAE(m/s)');
xlabel('Step');

%% close loopģʽ��ʵ��-��ֵ����,��ѵ�����е�Ч��,���ǶԲ��Լ��Ͳ�һ������,���Ըо����ǲ��ܼ���
% es_narx=[es(end,1:ri-1);es(:,2:ri)];
% qx_narx=[qxs(end,1:ri-1);qxs(:,2:ri)];
% 
% for i=1:96
%     for j=1:ri-1
%         phInputs_c=num2cell(qx_narx(i:i+1,j)');
%         if i==1
%             PhTargets_c=num2cell(es_narx(1:2,j)');
%         else
%             PhTargets_c=[yp1(i-1,j),num2cell(es_narx(i+1,j))];
%         end
%        [p1,Pi1,Ai1,t1(i,j)] = preparets(narx_net_closed,phInputs_c,{},PhTargets_c);
%        yp1(i,j) = narx_net_closed(p1,Pi1,Ai1);
%     end
%     error_narx2(i,:)=cell2mat(yp1(i,:))'-cell2mat(t1(i,:))';
%     m_narx(i,:)=mean(error_narx2(i,:));
%     yp1(i,:)=num2cell(cell2mat(yp1(i,:))-m_narx(i,:));
%     error_narxm(i,:)=error_narx2(i,:)-m_narx(i,:);
% end
% 
% for i=1:96
%     e_narxm(i,:)=mae(error_narxm(i,:));
% end
%% ������Ԥ��
load data80_test
count90=0;
count80=0;
count70=0;
for i=1:ri-1
    for j=1:96
        s90=interp1(qwb{1,j}(:,1),qwb{1,j}(:,2),qx_narx(j,i));
        x90=interp1(qwb{1,j}(:,1),qwb{1,j}(:,3),qx_narx(j,i));
        s80=interp1(qwb{1,j}(:,1),qwb{1,j}(:,4),qx_narx(j,i));
        x80=interp1(qwb{1,j}(:,1),qwb{1,j}(:,5),qx_narx(j,i));
        s70=interp1(qwb{1,j}(:,1),qwb{1,j}(:,6),qx_narx(j,i));
        x70=interp1(qwb{1,j}(:,1),qwb{1,j}(:,7),qx_narx(j,i));
        cha90(j,i)=qx_narx(j,i)+s90-max(qx_narx(j,i)+x90,0);
        cha80(j,i)=qx_narx(j,i)+s80-max(qx_narx(j,i)+x80,0);
        cha70(j,i)=qx_narx(j,i)+s70-max(qx_narx(j,i)+x70,0);
        if cfs(j,i+1)<=qx_narx(j,i)+s90 && cfs(j,i+1)>=qx_narx(j,i)+x90
            count90=count90+1;
        end   
        if cfs(j,i+1)<=qx_narx(j,i)+s80 && cfs(j,i+1)>=qx_narx(j,i)+x80
            count80=count80+1;
        end
        if cfs(j,i+1)<=qx_narx(j,i)+s70 && cfs(j,i+1)>=qx_narx(j,i)+x70
            count70=count70+1;
        end
    end
end
cha90(isnan(cha90))=0;
cha80(isnan(cha80))=0;
cha70(isnan(cha70))=0;
cha90m=mean(mean(cha90));
cha80m=mean(mean(cha80));
cha70m=mean(mean(cha70));

for i=1:ri-1
    mae_narx(i)=mae(error_narx(:,i));
end
[dayU,i]=sort(mae_narx);
i=i(cixu);
for j=1:96
    s90b(j,:)=interp1(qwb{1,j}(:,1),qwb{1,j}(:,2),qx_narx(j,i));
    x90b(j,:)=interp1(qwb{1,j}(:,1),qwb{1,j}(:,3),qx_narx(j,i));
    s80b(j,:)=interp1(qwb{1,j}(:,1),qwb{1,j}(:,4),qx_narx(j,i));
    x80b(j,:)=interp1(qwb{1,j}(:,1),qwb{1,j}(:,5),qx_narx(j,i));
    s70b(j,:)=interp1(qwb{1,j}(:,1),qwb{1,j}(:,6),qx_narx(j,i));
    x70b(j,:)=interp1(qwb{1,j}(:,1),qwb{1,j}(:,7),qx_narx(j,i));     
end
A=qx_narx(:,i)+x70b;
A(A<0)=0;
B=qx_narx(:,i)+x80b;
B(B<0)=0;
C=qx_narx(:,i)+x90b;
C(C<0)=0;
figure;
plot(qx_narx(:,i),'k','linewidth',1.2);
hold on
plot(qxs(:,i+1),'--r','linewidth',1.2);
hold on
plot(cfs(:,i+1),'m','linewidth',1.2);
hold on
ylabel('Wind Speed(m/s)');
a=[repmat(0.9,1,96),repmat(0.8,1,96)];
b=[repmat(0.8,1,96),repmat(0.7,1,96)];
hold on
fill([1:96 fliplr(1:96)],[(qx_narx(:,i)+s90b)' fliplr((qx_narx(:,i)+s80b)')],a,'EdgeColor','none');%��͸����,'facealpha',0.5
hold on
fill([1:96 fliplr(1:96)],[(qx_narx(:,i)+s80b)' fliplr((qx_narx(:,i)+s70b)')],b,'EdgeColor','none');
hold on
fill([1:96 fliplr(1:96)],[B' fliplr(A')],b,'EdgeColor','none');
hold on
fill([1:96 fliplr(1:96)],[C' fliplr(B')],a,'EdgeColor','none');
colorbar('SouthOutside')
% ylabel(colorbar,'p')
hold on
plot(qx_narx(:,i),'k','linewidth',1.2);
hold on
plot(qxs(:,i+1),'--r','linewidth',1.2);
hold on
plot(cfs(:,i+1),'m','linewidth',1.2);
hold on
ylim([0,12.5]);
ylabel('Wind Speed(m/s)');
set(gca,'FontName','Times New Roman','FontWeight','bold','Fontsize',24);
title('80m (Jun.11)','FontWeight','bold','Fontsize',24);
set(gcf,'position',[0,0,800 600]);
% legend('Corrected one-day-ahead WSF','NWP','Actual measured');
% save data80_test